In [4]:

    

import os
import knmi
import scipy
import netCDF4
import waterbase
import pandas as pd
import cPickle as pickle
from datetime import datetime

root = '/Users/hoonhout/GitHub/aeolis-models/sandmotor/'

t0 = datetime(2011,8,1)
t1 = datetime(2012,8,1)

key = '20110803'


    

In [2]:

    

# read bathy
with open('/Users/hoonhout/Checkouts/Zandmotor/morphology/JETSKI/bathy_gridded.pkl', 'rb') as fp:
    grid = pickle.load(fp)


    

In [5]:

    

def rotate(x, y, phi, origin=(72502.0, 452074.0)):

    x = np.asarray(x)
    y = np.asarray(y)
    phi = phi / 180. * np.pi

    R = np.asarray([[np.cos(phi), -np.sin(phi)],
                    [np.sin(phi),  np.cos(phi)]])

    XY = np.concatenate((x.reshape(-1,1) - origin[0],
                         y.reshape(-1,1) - origin[1]), axis=1)

    XY = np.dot(XY, R)
    xr = XY[:,0].reshape(x.shape)
    yr = XY[:,1].reshape(y.shape)
    
    return xr, yr


# select and rotate bathy
alpha = 48.
X, Y = rotate(grid['X'], grid['Y'], alpha)
Z = grid['Z'][key]


    

In [33]:

    

# interpolate bathy to computational grid and write to file

d = 50
x = np.arange(-500,  1000, d)
y = np.arange(-1500, 2500, d)

Xr, Yr = np.meshgrid(x, y)
np.savetxt(os.path.join(root, 'x.txt'), Xr)
np.savetxt(os.path.join(root, 'y.txt'), Yr)

for i, key in enumerate(sorted(grid['Z'].keys())):
    Z = grid['Z'][key]
    
    ix = ~np.isnan(Z.flatten())
    Zr = scipy.interpolate.griddata((Y.flatten()[ix],
                                     X.flatten()[ix]),
                                    Z.flatten()[ix],
                                    (Xr.flatten(),
                                     Yr.flatten()), method='nearest').reshape(Xr.shape)

    np.savetxt(os.path.join(root, 'z_%d.txt' % i), Zr)

print Zr.shape


    

    




(80, 30)

In [7]:

    

# plot computational bathy

fig, axs = subplots()
axs.pcolormesh(Yr, Xr, Zr, cmap='Spectral', vmin=-10, vmax=10)
axs.set_xlim((-1500, 2500))
axs.set_ylim((-500, 1000))
axs.set_aspect('equal')


    

In [11]:

    

# read wind time series
meteo, legend = knmi.read_uurgeg(os.path.join(root, 'uurgeg_330_2011-2020.txt'))
ix = (meteo['DD'] > 360.) | (meteo['FH'] >= 900)
meteo = meteo.ix[~ix]


    

In [51]:

    

# write wind time series
m = meteo[t0:t1].filter(['FH','DD'])
t = np.asarray([[x.total_seconds() for x in m.index - m.index[0]]]).T

m['FH'] = m['FH'].divide(10.)
m['DD'] = m['DD'].subtract(alpha)

np.savetxt(os.path.join(root, 'wind.txt'), np.concatenate((t, m.as_matrix()), axis=1))


    

In [14]:

    

# read and write water levels
wb = waterbase.Waterbase()
wl = wb.get_data('1', 'HOEKVHLD', t0, t1)
t = np.asarray([[(x - wl['datetime'][0]).total_seconds() for x in wl['datetime']]]).T
wl = np.concatenate((t, wl['value'].reshape((-1,1))/100.), axis=1)

np.savetxt(os.path.join(root, 'tide.txt'), wl)


    

In [18]:

    

# read and write wave heights
Hs = wb.get_data('22', 'Euro platform', t0, t1)
t = np.asarray([[(x - Hs['datetime'][0]).total_seconds() for x in Hs['datetime']]]).T
Hs = np.concatenate((t, Hs['value'].reshape((-1,1))/100.), axis=1)

np.savetxt(os.path.join(root, 'waves.txt'), Hs)


    

In [50]:

    

def create_plot(par, data):
    
    # remove inf's
    ix = np.isinf(data)
    data[ix] = np.nan
    
    # remove extra dimensions
    if data.ndim == 4:
        data = data[:,:,0,:] # pick top layer
    if data.ndim == 3:
        data = data.sum(axis=-1) # sum over fractions
    
    # determine data range
    mn = np.nanmin(data[1:-1,1:-1])
    mx = np.nanmax(data[1:-1,1:-1])
    mx = np.maximum(np.abs(mn), np.abs(mx))

    if mn > -1e-10:
        mn = 0.
        cmap = 'Reds'
    else:
        mn = -mx
        cmap = 'bwr'

    # create plot
    fig, axs = subplots(figsize=(10,4))
    p = axs.pcolormesh(y, x, data, cmap=cmap, vmin=mn, vmax=mx)
    axs.set_aspect('equal')
    axs.set_title(par)
    fig.colorbar(p)

    # save plot
    fig.savefig(os.path.join(root, 'plots', '%s.png' % par))

    
n = 10 # number of time steps
with netCDF4.Dataset(os.path.join(root, 'aeolis.nc'), 'r') as ds:
    t = ds.variables['time'][:] / 3600.
    x = ds.variables['x'][:,:]
    y = ds.variables['y'][:,:]
    f = ds.variables['fractions'][:]
    
    nx = x.shape[1]
    ny = y.shape[0]
    
    for t in np.linspace(1, len(t)-1, n):
        
        t = int(np.round(t))
    
        # bathy difference
        data = ds.variables['zb'][...]
        create_plot('dz_t=%d' % t, data[t,:,:] - data[0,:,:])

        # mean grain size
        data = ds.variables['mass'][...]
        gs = ds.groups['settings'].getncattr('grain_size') * 1e6
        create_plot('d50_t=%d' % t,
                    np.asarray([np.average(gs, weights=ds.variables['mass'][t,j,i,0,:])
                                for j in range(ny) for i in range(nx)]).reshape((ny, nx)))
        
        # other variables
        for par in ['Ct', 'Cu', 'pickup', 'zb', 'zs']:
            create_plot('%s_t=%d' % (par, t), ds.variables[par][t,...])


    

In [51]:

    




    

In [17]:

    




    

In [ ]: